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We present a method which computes many-electron energies and eigenfunctions by a full con- 
figuration interaction which uses a basis of atomistic tight-binding wave functions. This approach 
captures electron correlation as well as atomistic effects, and is well suited to solid state quantum 
dot systems containing few electrons, where valley physics and disorder contribute significantly to 
device behavior. Results are reported for a two-electron silicon double quantum dot as an example. 



I. INTRODUCTION 

Semiconductor devices have traditionally been simu- 
lated and understood using semiclassical band theory. 
An increasing number of today's nanoscalc semiconduc- 
tor devices, however, require more precise treatment of 
their quantum mechanical nature. Particularly within 
the area of quantum computation, capturing quantum 
effects is crucial to understanding the behavior and op- 
eration of devices. 

Simulation of quantum-mechanical semiconductor de- 
vices presents several challenges. The many degrees of 
freedom present due to the semiconductor lattice and the 
interactions between electrons quickly lead to intractable 
problems. Effective mass theory significantly reduces the 
number of degrees of freedom by approximating the ef- 
fects of the lattice with several effective parameters, but 
as a result is unable to capture atomistic effects such as 
interface roughness, lattice miscuts, defects, alloy disor- 
der, and valley splitting without introducing adjustable 
parameters. 

There is also a question of how approximately the 
electron-electron interactions can be treated in a quan- 
tum device. In the many-electron devices used in classical 
computing the interactions among electrons can be well 
approximated using Fermi liquid theory or by adding ef- 
fective potential terms, and in any case the precise wave 
function of the many-electron system is not nearly as 
important as the electron density. Devices proposed for 
quantum computation, however, often rely on manipu- 
lation of a wave function of few electrons and accurate 
modeling of electron correlation is essential. 

There are well-established methods to separately sat- 
isfy the demands to capture atomistic effects and 
electron-electron interactions. Tight binding methods^, 
which include orbitals for each atom of the lattice, are 
well-suited for modeling atomistic effects in crystalline 
materials. Although this method scales well in speed and 
memory with the device volume enabling several million 
atom systems to be solved rapidly^, it is essentially a 
single electron theory. Many-electron effects can only be 
captured in an approximate way by mean-field methods 
such as charge self-consistent iterations with the Poisson 
equation^ or a reduced basis Hartree Fock method^. Ex- 
act many-body exchange-correlations and wave functions 
are difficult to capture under this framework, especially 



for systems with a large number of atoms. 

The configuration interaction (CI) method^ is poten- 
tially one of the most precise ways of capturing electron- 
electron interactions. In it, the full quantum Hamiltonian 
is exactly diagonalized in a basis of multi-electron states. 
Since its computational requirements grow rapidly with 
particle number, the CI method is limited to systems 
with small numbers of electrons. In this paper, we de- 
scribe a method which integrates empirical tight bind- 
ing and full configuration interaction methods. We also 
present an algorithm which achieves a substantial speed 
up in the evaluation of the Coulomb and Exchange in- 
tegrals, which is the principle bottleneck of past at- 
tempts to integrate the two methods^—. The result 
is a tool which satisfies both of the issues discussed 
above simultaneously. In particular, it is capable of ac- 
curately modeling few-electron quantum devices involv- 
ing several million atoms while including atomistic ef- 
fects. Prior work has used tight-binding configuration 
interation methods for smaller systems such as molecules 
and nano crystals^ - — but employs different algorithms 
which are less well suited to the large number of atoms 
that must be considered quantum dot devices. 

As this method is motivated by the need to simu- 
late quantum dot systems to be used as quantum bits 
(qubits), we study as an example a silicon double quan- 
tum dot (DQD) qubit which encodes quantum informa- 
tion in the lowest energy two-electron singlet and triplet 
state a 16 ' 17 . Understanding and quantitatively estimating 
effects arising from the multi- valley nature of silicon, as 
well as the disorder which occurs at heterostructure in- 
terfaces, is crucial to an accurate understanding of the 
device operation. Furthermore, because quantum infor- 
mation is stored in two-electron states, precise account- 
ing of the interaction between these electrons is required. 
By design the simulation goals for this system align with 
our method's capabilities. We note that more approxi- 
mate methods, such as effective mass theory, may con- 
sider the effects which arise from the lattice but must 
then also include adjustable parameters which describe 
these effect o 18 ' 19 . In our method, such adjustable param- 
eters are not necessary because atomistic quantities can 
be computed directly. 

In the sections that follow, we describe the method and 
then present results for our example system. 
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II. METHOD 

A. Overview 

Our approach begins with an empirical tight binding 
solver, NEM03D2i2£ (Nanoelectronic Modeling Tool), 
which models a system as a collection of atoms, each 
with a specified number of orbitals, within the tight 
binding approximation. This accounts for atomic prop- 
erties of the material being simulated, such as valley 
coupling, defects, lattice strain, and interface roughness. 
The single-particle wave functions and energies gener- 
ated by NEM03D are input to a configuration interac- 
tion method which diagonalizes the many-body Hamilto- 
nian in the basis of Slater determinants constructed from 
NEM03D wave functions. Some of the main features of 
the two methods arc described below. 



1. Empirical tight binding method in NEM03D 

In the empirical tight-binding method, the Hamilto- 
nian of the lattice is described in a basis of atomic orbitals 
localized in each atom. As a result, the electronic wave 
function is expressed as a linear combination of these 
atomic orbitals^. Following Slater and Kostcr— , best re- 
sults are obtained if the Hamiltonian matrix elements 
are optimized to fit critical features of the bulk band 
structure, such as band gaps, band minima locations, 
band symmetries, effective masses, and so on. Once a set 
of tight-binding Hamiltonian parameters are found for 
a host unit cell, they can used to compute the detailed 
electronic structure of any superlattice of the material 
along with an added external potential. This is a stan- 
dard technique referred to as the empirical tight binding 
method. 

In this work, we utilize the more complete 10 band 
sp 3 d 5 s* nearest neighbor model for silicon to represent 
the Hamiltonian^. The model parameters are optimized 
by an advanced genetic algorithm procedure^, and are 
well established in literature. The method, in general, 
is able to capture the detailed electronic structure of 
systems of different dimensions, such as bulk, quantum 
wells, wires and dots with respect to experiments^— . 
The approach outlined in this paper is not specific to 
this band model, but can be applied generally to other 
band models and other materials. If the corresponding 
20 band spin model is used, spin-orbit and Zeeman split- 
tings can be captured directly from the eigen spectrum. 

Bandstructure: Since the tight binding method is a full 
band-structure method, quantum states near all conduc- 
tion or valence band valleys can be captured. Effects 
such as valley-splitting can be obtained without any ad- 
ditional adjustable parameters. 

Strain: If the device has at least two materials of dif- 
ferent lattice constants, the relaxed structure is obtained 
by minimizing the strain energy by the valence force field 
method of Keating^. This captures the effects of inho- 



mogencous strain distribution near hetero-structure in- 
terfaces, where quantum dot wave functions usually re- 
side. Alternately, homogeneous strain can be simulated 
by setting the lattice constants to the relaxed values 
given by analytic expressions for various alloys^. The 
strain induced modifications to the Hamiltonian matrix 
elements have been worked out in^^. 

Alloys: Alloyed materials such as Sii-^Ge^ can be 
treated both atomistically 3 ^ or from a virtual crystal ap- 
proximation. An atomistic treatment allows the investi- 
gation of realistic systems relevant to experiments where 
disorder causes sample-to-sample variations in the mea- 
sured electronic properties. 

External fields: Electrostatic potential either gener- 
ated analytically or in TCAD (Technology Computer 
Aided Design) tools can be interpolated on the atomistic 
lattice, and can be included in the Hamiltonian. Hence, 
it is possible to study electronic states as a function of 
gate bias, which is the basis of all quantum nanoelectron- 
ics. Magnetic fields can be included both in symmetric 
or asymmetric gauges depending on the system. 

Geometry and interfaces: Complex device geometries 
can be described as a collection of simple shapes such as 
cuboids, cylinders, spheres, and so on. In general, re- 
alistic device geometries described in other tools can be 
imported, as long as the materials in the simulation do- 
main are either crystalline or are described by a virtual 
crystal model. Hetero-structure interfaces are modeled 
naturally as the atom positions are defined. Hydrogen 
passivated interfaces have also been modeled 3 -^. More- 
over, miscut interfaces and surface roughness can also be 
defined atomistically, and used in the simulations^. 

Solution: The full Hamiltonian is solved by a paral- 
lel Lanczos algorithm to obtain the desired number of 
cigenstates in a specified energy range of interest. If de- 
generate or nearly degenerate eigen states are required, 
a Block Lanczos algorithm is used with a block size cho- 
sen to meet the degree of resolution of the eigen states 
required 3 ^. 

Capability: NEM03D is capable of solving large device 
volumes. Single electron wave functions of a few coupled 
quantum dots can be solved fairly quickly. The code 
has been demonstrated so far to solve 52 million atoms 
(101 nm 3 ) for electronic structure and over 100 million 
atoms for strain relaxation with very impressive scaling 
in memory and solution time^S. 



2. Configuration Interaction 

In the CI method, a basis of many electronic config- 
urations represented by Slater Determinants^ are con- 
structed using the single electron wave functions. Each 
configuration corresponds to electrons occupying partic- 
ular single electron energy levels. A Slater Determinant 
in which an electron has been promoted to an excited 
orbital from the ground state represents a singly excited 
configuration of the system. For an N electron system, 
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the highest excitation of the system consists of all ./V elec- 
trons occupying a subset of the excited orbitals. If the 
basis set for the N electron CI consists of all excitations, 
it is termed a full CI (FCI). In cases where a full CI is 
intractable, the number of excitations considered can be 
limited. 

The iV-electron many body Hamiltonian can be 
separated into single- and two-particle parts, % = 
J2iLi'H-i P + J2i<jUif- 1i\ p is a diagonal matrix in 
the basis of the single electron states ipi, with eigen en- 
ergies Ei, which are obtained from NEM03D. The two- 
particle part is simply the Coulomb interaction, Hif = 
e 2 / 1 ?i — Tj\. This Hamiltonian is evaluated in the basis 
of the Slater Determinants to obtain the CI matrix. The 
technical details of an FCI method based on Gaussian 
orbitals has been described in our earlier work^. The 
CI matrix is diagonalized using a Lanczos sparse matrix 
algorithm. When the Hamiltonian commutes with total 
spin and/or the total of a component of the spin {e.g. to- 
tal spin ^-component), the CI matrix is first block di- 
agonalized according to these symmetries and then each 
block is diagonalized separately. 



B. Procedure 

The following is a step-by-step overview of the method. 

1. System/device is specified within NEM03D using 
regions of different materials, complex geometries, 
and external potential. Domains for the computa- 
tion of strain and electronic structure (not neces- 
sarily the same) are chosen. 

2. NEM03D computes the n lowest single particle en- 
ergies Ei and wave functions ipi of the system, pos- 
sibly performing lattice relaxation first to account 
for strain. 

3. CI module discretizes each ipi on a sub-atomic 3D 
rectangular grid by assuming analytic Slater-type 
orbital forms for the unknown basis functions used 
by NEM03D. This step is essentially performs a 
low-pass filter on the wave function. 

4. The Fourier transform is taken of every possi- 
ble product of two discretized single particle wave 
functions, resulting in n(n + l)/2 "pair products" 
Hjj(r) = FT (ipi(r)ipj(r)) (discretized on the rect- 
angular fc-spacc lattice reciprocal to the spatial 
grid). It is possible, at this stage, to only save por- 
tions of each pair product function with large mag- 
nitude (i.e. large spectral weight, determined by a 
cutoff), just as sparse matrices store only nonzero 
elements. This sparse storage in the Fourier do- 
main is similar to JPEG image compression, and 
can significantly reduce the storage requirements 
for the method. 



5. The CI Hamiltonian is constructed using the single 
electron energies Ei obtained from NEM03D and 
the two-particle Coulomb interaction terms ~H^f = 
e 2 /\fi — fj \ evaluated using the pre-computcd pair 
products as described in the detailed discussion be- 
low. . 

6. CI module diagonalizes the Hamiltonian matrix in 
the basis of Slater determinants of the ipi, and 
outputs the many-electron energies and wave func- 
tions. 



C. Approximations 

Several approximations/assumptions are used in this 
method. The first is the choice of the atomic orbitals for 
the tight-binding wave functions. As we discussed ear- 
lier, in the empirical tight-binding method, the Hamil- 
tonian matrix elements comprising on-site energies and 
nearest-neighbor interactions are optimized numerically 
to fit the bulk band structure of the host material. The 
symmetries of these orbitals determine which terms in 
the Hamiltonian are non- vanishing. No explicit forms of 
these orbitals are assumed. 

However, in order to evaluate the Coulomb and Ex- 
change integrals within the CI, the spatial forms of these 
orbitals need to be specified. There is no other way to 
do this than to arbitrarily choose such a set of orbitals, 
as long as they adhere to the correct orbital symmetries 
and represent the atomic orbitals of the host to a certain 
degree. These basis sets cannot be related to the op- 
timized tight-binding matrix elements in any way. The 
choice of the basis set can only be justified by comparing 
the many electron levels with those computed with other 
basis sets, and studying the sensitivity of the results on 
the choice. Benchmarking with experimental data can 
be another justification for the choice. 

A prior work had developed a real space CI method 
based on tight-binding wave functions. Since the 
Coulomb and Exchange integrals were evaluated in real 
space rather than in momentum space, as we do here, 
the method is limited to dots of only 3-5 nm diame- 
ter. However, the work investigated the sensitivity of 
the many-electron states for a number of basis sets, such 
as Slater Type Orbitals, orthogonal and non-orthogonal 
Gaussian Type Orbitals. The work also showed that as 
the dot radii increased, the Coulomb and Exchange inte- 
grals became less sensitive to the choice of the basis set. 
In fact, beyond dots of radii 1.5-2.0 nm, the integrals 
were found quite reliable and representative of the actual 
values. This can be understood from the form of the in- 
tegrals. If the Coulomb interaction is between two points 
a large distance apart, then the atomic orbitals appear as 
point-like charges, and their shapes do not influence the 
computation. Fortunately, the quantum dot wave func- 
tions we are considering span about 30 nm in each of the 
lateral dimensions. For such systems, the choice of the 
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atomic basis set can be completely arbitrary. We have 
chosen STOs as they can be obtained easily for each ma- 
terial from the Slater rules^, and are representative of 
the atomic scale screening of the materials. 

We choose the following Slater-type forms for the two 
lowest s-, three lowest p-, and five lowest d-type symme- 
try orbitals, since they are good approximations to the 
orbitals of an isolated atom. Expressions for these or- 
bitals for silicon are given in Table |TJ The need to choose 
a set of basis functions will arise whenever an empiri- 
cal tight binding and configuration interaction techniques 
are coupled, and in this sense, the choice of a certain 
atomic basis is a necessary assumption for this method. 



Orbital type 


wave function 


s 


1.30171r 2 e- 1 - 38r 


P 


1.30171r 2 e- 1 ' 38r yi m (6l,V) 


d 


0.037268r 2 e- a5r y 2 m(6',^) 


s* 


0.00337r 27 e- - 39 '' 



TABLE I: Basis functions for silicon, as given by the Slater 
rules,— where Y; m is a spherical harmonic and r is assumed 
to be in units of Bohr radii.— 

A second approximation is the discretization of the 
tight-binding wave function, which is a continuous func- 
tion (but dependent on the specific choice of atomic basis 
functions), on a discrete uniform mesh (grid). The con- 
stant spacing of this mesh determines a high-frequency 
cutoff for the resulting discrete function. In order to cap- 
ture atom-scale detail in the wave function, this cutoff 
must be large enough to capture the highest frequencies 
with appreciable weight of the Slater-type basis functions 
defined above. This roughly corresponds to there being 
many mesh points on the scale of the basis functions, 
so that the shape of each basis function is well approxi- 
mated using the given mesh. When this criteria is met, 
the discretization represents only a minor approximation, 
especially since the exact form of the continuous basis 
functions are not known. The reason for this discretiza- 
tion is performance-motivated, as will be described be- 
low, and other ways of coupling TB and CI methods do 
not require it. 

The remaining approximations are those intrinsic to 
any CI method: a subset of the single particle functions, 
in our case the n with lowest energy, are used to form 



the many-electron basis states. In the appendix we ana- 
lyze the convergence of the CI results with n in a double 
quantum dot device. The number of many-electron ba- 
sis states can be reduced by limiting the solution space 
to n ex excitations relative to a given reference state (see 
Ref£ for details on CI methods), be between 1 

and the number of particles N, at which point the CI is 
termed a "full CI" . In the example shown in section IIIII 
below, n = N = 2 so the CI is a full CI. 



D. Performance 

The CI portion of the method, whether full or trun- 
cated, scales in the usual way with the number of parti- 
cles N and size of the single particle basis n. In the full-CI 
case, this scaling is exponential in N and n. This is due to 
the size of the Hamiltonian matrix, whose construction 
and partial diagonalization usually dominate the com- 
putation. When the single particle basis functions are 
defined on atomistic length scales, however, there is an- 
other potential bottleneck in the calculation: the eval- 
uation of the two-electron interaction term e 2 /|fi — r%\. 
In cases where there are many atomic sites but few elec- 
trons, the computational resources required to construct 
the matrix H 2P (cf. step [5] above) can dominate even the 
diagonalization of the CI matrix. Low-electron quantum- 
dot devices are a prime example of systems which have 
few electrons in regions containing many atomic sites. 

It is to reduce the cost of evaluating the elements of 
V. 2P that we perform the discretization in step [3] and 
compute Fourier transforms of "pair products" , ILj , in 
step 21 Let us denote the number of atomic sites as M. 
Since the number of mesh points at which we evaluate 
each single-particle wave function is proportional to the 
number of atomic sites, and since only a constant num- 
ber of neighboring sites are needed to evaluate a wave 
function at a point (due to the localized nature of the 
atomic basis), the cost of evaluating each wave function 
on the discrete mesh is proportional to M, Multiplying 
two wave functions together and taking the fast Fourier 
transform of the result requires time 0(M ln(M)), so that 
in total the time required to compute the n(n + l)/2 pair 
products is 0(nM + n 2 Mln(M)) = 0(n 2 M\n(M)). Af- 
ter all Hij are computed, the matrix elements of T~L 2P are 
computed as 
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\fi - r 2 \ 



ipairDipbirli)-, 



\ri - r 2 \ 



PaM)-. 



\n - r 2 \ 

4vT, 



fj e iA 3 (n-r- 2 )^ 3 
k 3 



(2ir) h / U ac (-k)—n M (k)dk. 



r 



(i) 



(2) 



n bd (fc 2 )e lfc2r2 dfc 2 ) dr ± dr 2 (3) 

(4) 



In the final line, we have an expression for an element 
of the H 2P matrix, in the basis of products ipiipj, as a 
single (3D) integral over reciprocal space, with functions 
IT ac and Tlbd prccomputed as described. As there are 
the same number of reciprocal space points as those of 
the real space mesh, the time to perform each of the 
0(n 4 ) integrals is proportional to M, and thus the time 
to construct the T-L 2P matrix is 0{n A M). 

This can be compared with the more standard ap- 
proach whereby each element {ip a ipb\ r^z^i IV'cV'd) is com- 
puted as a double integral (or double sum over prod- 
ucts of basis function coefficients) which requires time 
0(n 4 M 2 ) to construct the H 2P matrix. This extra factor 
of M, which can be of order 10 6 — 10 7 , makes this direct 
method much slower than the one we describe (note the 
0(n 2 M ln(M)) time required to pre-compute the Iljj is 
insignificant). 



III. EXAMPLE: TWO-ELECTRON DOUBLE 
QUANTUM DOT 

We now apply the method to compute the two-electron 
states of a silicon DQD, an important implementation of 
a solid-state qubit. It was shown that the lowest single- 
triplet (i.e. unpolarized triplet To) of a DQD form a two- 
level system ideal to encode quantum informatio n 16 ' . 
Recent experiments^ have demonstrated spin initializa- 
tion, coherent manipulation, and readout of this two- 
level system in a GaAs DQD. Due to the intrinsic nu- 
clear spin associated with the host atoms in GaAs, the 
spins of the dot bound electrons decohere faster, even on 
the scale of nanoseconds to microseconds as t 2 * measure- 
ments suggests^ 6 -. Since silicon has a very low concen- 
tration of isotopes with nuclear spin and can be enriched 
to obtain high purity silicon, quantum confined electrons 
in silicon can have very long spin coherence times - even 
of order seconds 3 ^. Silicon DQDs are therefore an ideal 
choice for long-lived qubits. However, silicon systems are 
more prone to disorder at the Si-Si02 interface arising 
from charge defects or surface roughness. In addition, the 
conduction band valley degeneracy adds an additional de- 
gree of complexity in the electronic structur o 18 ' 26 ' 37 " — . 

Our method is ideal to study a silicon DQD in the pres- 



ence of such imperfections. We show 2-clcctron energy 
level structure which includes valley effects and agrees 
with effective mass calculations. The valley splitting ap- 
pears without introducing an explicit valley term in the 
Hamiltonian, and we see that this splitting varies as ex- 
pected with the vertical electric field strength. 

Consider a two-electron double quantum dot given by 
the Hamiltonian 



(Pi - eAf 



i=l 



2m 



+ V(f I ) + g^ B S I -B + 



i<j 



(5) 

where r*j and pi are the position and momentum, respec- 
tively, of the i th electron, V is the electrostatic potential, 
to is the (bare) electron mass, and k is the silicon di- 
electric constant. A vector potential A determines the 
magnetic field B = V x A, which we set to zero for the 
cases considered here. The potential is the sum of the 
potential due to the atomic lattice and an external po- 
tential, V = Viatt + V ex t- We idealize the DQD potential 
Vext as the minimum of two parabolic dots with a vertical 
electric field F z , 

Vext(x, y) = a [min ((x - L) 2 + e, (x + L) 2 ) + y 2 ] +F z z . 

(6) 

The parameters e, L, and a, correspond to the bias, inter- 
dot distance, and a measure of the confining well poten- 
tial or dot size, respectively. Vi at t is set by the mate- 
rials); in this case silicon with a Hydrogcn-passivatcd 
surface. 

Qubit manipulation in the DQD system relies on a 
voltage controlled exchange splitting between the singlet 
and the triplet state. For this reason, one dot is detuned 
relative to the other by a gate bias e such that the elec- 
trons in the dots make an adiabatic transition from a 
(1,1) occupation to a (0,2) occupation. A (1,1) occupa- 
tion refers to one electron in each dot, a system in which 
the overlap between the electronic wave functions at the 
tunnel barrier between the dots is small or negligible, re- 
sulting in a very small exchange energy. As the dots are 
swept into the (0,2) occupation, which refers to both elec- 
trons in the same dot, the system undergoes a transition 
to a high overlap system with increased exchange energy. 
As the detuning is performed adiabatically, regions of the 
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FIG. 1: Energy levels of a silicon (passivated-surface) DQD 
as a function of detuning parameter e. L = 30 nm, a = 
0.0001 eV/nm 2 , and F z = 5mV/nm. 
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FIG. 2: Energy levels of a silicon (passivated-surface) DQD 
as a function of detuning parameter e. L = 30 nm, a = 
0.0001 eV/nm 2 , and F z = 20mV/nm. 



exchange curve in between these two extremums are also 
accessible. 

While this control sequence is well-established in a 
GaAs DQDil, a Si DQD has additional valley degrees 
of freedom the impact of which depends on atomic scale 
effects such as miscuts, surface roughness, alloy disorder, 
and applied fields. The electronic states in a Si DQD 
therefore appears quite different from those in the GaAs 
system. We now present results for the DQD system 
which include these effects to illustrate the capabilities 
of this method. 



A. Valley splitting 

The splitting in energy between different valleys in 
semiconductor arises from a material's atomic properties 
and crystal lattice symmetries. In order to compute val- 
ley splitting from first principles, one must have an atom- 
istic model. Figure Q] shows the unpolarized (S z = 0) 
two-electron energy levels of a DQD parametrized by 
L = 20 nm, a = 0.0001 eV/ nm 2 , F z = 5mV/nm. We 
observe that the energy levels look like those of a DQD 
in a single valley material, but with a group of four states 
for each one state of the single valley case. This is due 
to the extra valley degree of freedom, as explained in 
RefJ£. The separation between the central line and ei- 
ther the upper or lower line is the "valley splitting" , and 
the spacings of the avoided crossings give a measure of 
the intra- and inter- valley coupling. We emphasize that 
these quantities are obtained without adjustable valley 
parameters, in contrast to effective mass methods. 

Figure [5] shows how the energy levels of Fig. [1] 
change upon increasing the vertical electric field F z to 
20 mV/ nm. In particular, the valley splitting increases 



with larger F z as expected. 

B. Atomic perturbations 

Perturbations at atomistic length scales are nearly 
ubiquitous in solid state systems. In this section we con- 
sider two types of atomistic perturbations, which we re- 
fer to as tilt and roughness. By tilt, we mean that the 
silicon-oxide interface is not exactly parallel to a crys- 
tallographic plane. This occurs due to miscuts in the 
silicon wafer and the miscut angle is typically one half 
to a few degrees. By roughness, we mean that the step 
edges at the silicon interface due to tilt are not straight 
but are located at a randomly varying positions. Since 
the QD electrons are confined at the interface between 
silicon and the barrier material, accurate models of in- 
terfaces are extremely important in electronic structure 
simulations. The TB method offers a natural way of deal- 
ing with tilt and surface roughness due to its atomistic 
nature. 

Miscuts are always present in wafers because it is usu- 
ally not possible to cut a wafer exactly along one of the 
crystallographic planes. Due to the miscut, the (001) sur- 
face of silicon reconstructs to form steps of monoatomic 
height and varying lengths. High precision scanning 
tunnelling microscope (STM) images have revealed the 
specific nature of these steps^ and models of the re- 
constructed surfaces have been developed from equilib- 
rium statistical mechanics based on the images4£ In this 
work, we have employed interface models of Refs. [26ll43l 
to demonstrate the advantages of atomistic description 
of the interfaces for qubit simulations. 

In particular, the (001) surface of silicon reconstructs 
to form rows of dimerized atoms. As a result, two types 
of steps are noticed in the STM images, one set of steps 
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(a) (b) 

FIG. 3: A depiction of the interface roughness present at a 
silicon (001) surface, which consists of alternating straight 
and rough steps due to dimerized atoms (see text). In (a) the 
real-space 3D profile of the surface is shown, and in (b) this 
view is projected onto the x-y plane. Such realistic disorder 
can be modeled atomistically without adjustable parameters 
by the technique described. 



running parallel and another perpendicular to the dimcr 
rows. The parallel steps are straight, whereas the per- 
pendicular steps have a high kink density, resulting in 
alternating straight and rough steps 42 A pictorial de- 
scription of these steps is shown in Fig. [31 with a 3D 
plot in (a) and a 2D plot in (b). Based on this rough- 
ness profile, the interface is generated in TB by adding 
or removing sets of the surface atoms ^ The addition of 
the CI technique now helps to understand the effect of 
these step disorder on the two-electron coupling essential 
in qubits. 

When the silicon interface is ideally tilted with no 
roughness, the valley splitting decreases relative to 
the smooth interface case due to interference between 
Bloch components of the wave function with different 
phase a 26 ' 38 . This effect is clearly seen by comparing 
Figs. [1] and |4j the latter of which has a rough inter- 
face tilted by 1 degree. We have chosen L = 20 nm in 
Fig. @] to better display the levels at low epsilon. The 
case where L = 30 nm, as in Fig. [T] shows a similar re- 
duction in valley splitting. The addition of roughness 
results in washing away some of the interference effects 
produced by ideal steps due to the randomness. This 
increases valley splitting compared to the ideally tilted 
case (without roughness). However, the valley splitting 
is still less than that of a flat surface with no roughness. 

IV. SUMMARY 

We have described a technique which couples tight 
binding and configuration interaction methods to make 
possible the accurate calculation of many-electron sys- 
tems in the presence of atomistic effects such as valley 
splitting or disorder. To illustrate the method's capabili- 
ties we apply it to a silicon double quantum dot contain- 
ing two electrons, where we find that varying the vertical 
electric field and adding interface tilt and roughness have 
the expected effects on the two electron states. 

This work was supported by the Laboratory Directed 
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FIG. 4: Energy levels of a silicon (passivated-surface) DQD 
as a function of detuning parameter e. The surface is tilted 
with respect to the crystallographic plane by 1 degree, and 
includes step roughness (the step edges due to the tilt are 
not straight). These atomistic properties reduce the valley 
splitting considerably compared to the smooth surface case 
(Fig. [2|. The slope of the lines relative to the smooth case 
is due to the vertical electric field not being perpendicular 
to the silicon surface and thereby contributing to the overall 
detuning of the DQDs. L = 30 nm, a — 0.0001 eV/ nm 2 , and 
F z = 20mV/nm. 
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Appendix A: Convergence 

As explained in the text, the configuration interaction 
(CI) algorithm takes as input a fixed number of single 
particle levels, n from which to generate Slater Deter- 
minant multi-electron basis functions. It is important 
to check convergence with respect to n in order to en- 
sure that the basis being used by the CI is adequate at 
approximating the full quantum Hilbert space of wave 
functions. 

Figures [5] and [5] show the energy levels and exchange 
(singlet-triplet) splitting, respectively, of a DQD for n = 
10, 14, and 18. We find here and in general for quantum 
dots with radii from 15 — 25 nm that results which use 
n = 10 arc sufficiently converged. The exchange energy 
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FIG. 5: Convergence of two-particle energies in a silicon DQD, 
given by the TB-CI method with respect to n, the number of 
NEM03D wave functions used by the FCI. DQD parameters 
L = 30 nm, a — 0.0001 eV/ nm 2 , and F z = hmV/mn. 



in Fig. [5] is taken to be the difference between the lowest 
singlet and unpolarizcd triplet states with the same val- 
ley character. (This means that at large values of e we 
consider the lowest energy singlet and the lowest triplet 
of the group of states with positive slope.) We expect 



when strong atomistic defect potentials are introduced 
to a quantum-dot system that a greater number of single 
particle levels will be required to obtain convergence. 
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FIG. 6: Convergence of the exchange energy in a silicon DQD, 
given by the TB-CI method with respect to n, the number of 
NEM03D wave functions used by the FCI. DQD parameters 
L — 30 nm, a = 0.0001 eV/ nm 2 , and F z = 5mV/nm. 
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